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Abstract 

We present a new model for the electricity spot price dynamics, which is able to cap- 
ture seasonality, low- frequency dynamics and the extreme spikes in the market. Instead 
of the usual purely deterministic trend we introduce a non-stationary independent in- 
crements process for the low-frequency dynamics, and model the large fluctuations by 
a non-Gaussian stable CARMA process. The model allows for analytic futures prices, 
and we apply these to model and estimate the whole market consistently. Besides 
standard parameter estimation, an estimation procedure is suggested, where we fit the 
non-stationary trend using futures data with long time until delivery, and a robust 
L 1 -filter to find the states of the CARMA process. The procedure also involves the 
empirical and theoretical risk premiums which - as a by-product - are also estimated. 
We apply this procedure to data from the German electricity exchange EEX, where 
we split the empirical analysis into base load and peak load prices. We find an overall 
negative risk premium for the base load futures contracts, except for contracts close to 
delivery, where a small positive risk premium is detected. The peak load contracts, on 
the other hand, show a clear positive risk premium, when they are close to delivery, 
while the contracts in the longer end also have a negative premium. 
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1 Introduction 



In the last decades the power markets have been liberalized world-wide, and there is a 
large interest for modelling power spot prices and derivatives. Electricty spot prices are 
known to be seasonally varying and mean-reverting. Moreover, a distinctive characteristic 
of spot prices is the large spikes that occur due to sudden imbalances in supply and demand, 
for example, when a large production utility experiences a black-out or temperatures are 
suddenly dropping. Typically in these markets, different production technologies have big 
variations in costs, leading to a very steep supply curve. Another characteristic of electricity 
is the lack of efficient storage possibilities. Many spot price models have been suggested 
for electricity, and we refer to Eydeland and Wolyniec [21] and Benth, Saltyte Benth and 
Koekebakker [3J for a discussion on various models and other aspects of modelling of energy 
markets. 

In this paper we propose a two-factor aritmethic spot price model with seasonality, which 
is analytically feasible for pricing electricity forward and futures contracts. The spot price 
model consists of a continuous-time autoregressive moving average factor driven by a stable 
Levy process for modelling the stationary short-term variations, and a non-stationary long- 
term factor given by a Levy process. We derive futures prices under a given pricing measure, 
and propose to fit the spot model by a novel optimization algorithm using spot and futures 
price data simultaneously. We apply our model and estimation procedure on price data 
observed at the German electricity exchange EEX. 

In a seminal paper by Schwartz and Smith |35j a two-factor model for commodity spot 
prices is proposed. Their idea is to model the short-term logarithmic spot price variations 
as a mean-reverting Ornstein-Uhlenbeck process driven by a Brownian motion, reflecting 
the drive in prices towards its equilibrium level due to changes in supply and demand. 
But, as argued by Schwartz and Smith [35J, there may be significant uncertainty in the 
equilibrium level caused by inflation, technological innovations, scarceness of fuel resources 
like gas and coal etc.. To account for such long-term randomness in prices, Schwartz and 



2 



Smith [35] include a second non-stationary factor being a drifted Brownian motion, possibly 
correlated with the short-term variations. They apply their model to crude oil futures traded 
at NYMEX, where the non-stationary part is estimated from futures prices, which are far 
from delivery. Mean-reversion will kill off the short-term effects from the spot on such futures, 
and they can thus be applied to filter out the non- stationary factor of the spot prices. 

This two-factor model is applied to electricity prices by Lucia and Schwartz [2E]- Among 
other models, they fit an arithmetic two-factor model with deterministic seasonality to elec- 
tricity spot prices collected from the Nordic electricity exchange NordPool. Using forward 
and futures prices they fit the model, where the distance between theoretical and observed 
prices are minimized in a least squares sense. 

A major critiscism against the two-factor model considered in Lucia and Schwartz [28] is 
the failure to capture spikes in the power spot dynamics. By using Brownian motion driven 
factors, one cannot explain the sudden large price spikes frequently observed in spot data. 
Multi-factor models, where one or more factors are driven by jump processes, may mend 
this. For example, Benth, Kallsen and Meyer-Brandis [2] suggest an arithmetic spot model 
where normal and spike variations in the prices are separated into different factors driven by 
pure-jump processes. In this way one may model large price increases followed by fast speed 
of mean- reversion together with a "base component" -behaviour, where price fluctuations are 
more slowly varying around a mean level. Such multi-factor models allow for analytic pricing 
of forward and futures contracts. 

A very attractive alternative to these multi-factor models are given by the class of continuous- 
time autoregressive moving-average processes, also called CARMA processes. These pro- 
cesses incorporate in an efficient way memory effects and mean-reversion, and general- 
ize Ornstein-Uhlenbeck processes in a natural way (see Brockwell [11J). As it turns out, 
(C)ARMA processes fit power spot prices extremely well, as demonstrated by Bernhardt, 
Kliippelberg and Meyer-Brandis [7J and Garcia, Kluppelberg and Miiller [23]. In [7J an 
ARMA process with stable innovations, and in [23] a CARMA(2,l)-model driven by a sta- 
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ble Levy process are suggested and empirically studied on power spot price data collected 
from the Singapore and German EEX markets, respectively. A CARMA(2,1) process may 
be viewed on a discrete-time scale as an autoregressive process of order 2, with a moving 
average order 1. By invoking a stable Levy process to drive the CARMA model, one is 
blending spikes and small innovations in prices into one process. We remark in passing that 
a CARMA dynamics has been applied to model crude oil prices at NYMEX by Paschke and 
Prokopczuk [30J and interest-rates by Zakamouline, Benth and Koekebakker |39j . 

We propose a generalization of the stable CARMA model of Garcia et al. [23] by including a 
long-term non-stationary factor being a general Levy process. The model allows for analytical 
pricing of electricity futures, based on pricing measures, which preserve the Levy property 
of the driving processes. More precisely, we apply an Esscher measure transform to the non- 
stationary part, and a transform which maps the stationary stable process into a tempered 
stable. Due to the semi-affme structure of the model, the futures and forward prices becomes 
explicitly dependent on the states of the CARMA model and the non-stationary factor. 

The CARMA-based factor in the spot model accounts for the short-term variations in prices 
and will be chosen stationary. By a CARMA model with a higher order autoregressive part 
we may include different mean-reversion speeds, such that we can mimic the behaviour of a 
multi-factor model accounting for spikes and base variations separately. The moving average 
part is necessary to model the observed dependence structure. The stable Levy process may 
have very big jumps, which then can explain spike behaviour in the prices. The smaller 
variations of the stable Levy process model the base signal in power prices. As it turns out 
from our empirical investigations using market price data from the EEX, the non-stationary 
long-term behaviour may accurately be modelled using a normal inverse Gaussian (NIG) 
Levy process. We filter out the non-stationary part from observing futures prices, which 
are far from delivery. The influence of the stationary CARMA factor is then not present, 
and the data shows a significant non-normal behaviour. This is in contrast to the choice 
suggested by Schwartz and Smith [35J and applied by Lucia and Schwartz [28J. Moreover, 
we find that a CARMA(2,1) model is accurately explaining the mean-reversion and memory 
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effects in the spot data. 

A novelty of our paper apart from the generalizing existing one and two factor models, is our 
estimation procedure. Lucia and Schwartz [25] propose an iterative algorithm for estimating 
their two-factor model to NordPool electricity data, where they minimize the least-squares 
distance between the theoretical and observed forward and futures prices to find the risk- 
neutral parameters. In order to find the theoretical prices, they must have the states of 
the two factors in the spot model accessible. Since these are not directly observable, they 
choose an iterative scheme, where they start with a guess on the parameter values, find 
the states minimizing distance, update parameters by estimation, find the states minimizing 
the distance etc. until convergence is reached. We propose a different approach, utilizing 
the idea in Schwartz and Smith [35] that the non-stationary factor is directly observable, 
at least approximately, from forward prices, which are far from delivery. We apply this to 
filter out the non-stationary factor. The CARMA-part is then observable from the spot 
prices, where seasonality and the non-stationary term is subtracted. Since we work with 
stable processes, which do not have finite second moments, L 2 -filters can not be used to 
find the states of the CARMA-process. We propose a simple L 1 -filter being more robust 
with respect to spikes in spot data to do this. The problem we are facing is to determine, 
what contracts to use for filtering out the non-stationary part. To find an optimal "time- 
to-maturity" which is sufficiently far from delivery, so that the futures prices behave as the 
non-stationary factor and at the same time provide a sufficiently rich set of data, we use an 
optimization algorithm, which minimizes the least square distance between the empirical and 
theoretical risk premium. In order to find the risk premia, we must have all the parameters 
of the model available, which in turn can only be found, if we know which futures contracts 
can be used for filtering the long-term factor. We implement an algorithm, which estimates 
all model parameters for futures contracts withs different times to maturities and minimizes 
the distance to the empricial risk premium. 

We apply our model and estimation scheme to data from the German EEX (European 
Energy Exchange) market where we use spot prices as well as futures prices of contracts 
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with a delivery period of one month. Our empirical studies cover both base load and peak 
load contracts, where base load contracts are settled against the average of all hourly spot 
prices in the delivery period. Peak load futures contracts are settled against the average of 
hourly spot price in peak periods of the delivery period. The peak load period is the period 
between 8 a.m. and 8 p.m., during every working day. As a first summary, we can say that 
the results for both base and peak load data are in general rather similar. However, the peak 
load data show a more extreme behaviour. The average risk premium decays when time to 
maturity increases, and is negative for contracts in the longer end of the futures curve. This 
points towards a futures market, where producers use the contracts for hedging and in return 
accept to pay a premium to insure their production, in accordance with the theory of normal 
backwardation. The risk premium is completely determined by the effect of the long-term 
factor, which induces a close to linear decay as a function of "time-to-maturity" . We see that 
for the base load contracts the risk premium in the short end of the curve is only slightly 
positive. The risk premium is negative for contracts starting to deliver in about two months 
or later. On the other hand, the peak load contracts have a clear positive risk premium, 
which turns to a negative one for contracts with delivery in about four months or later. The 
positive risk premium for contracts close to delivery tells us that the demand side (retailers 
and consumers) of the market is willing to pay a premium for locking in electricity prices as 
a hedge against spike risk (see Geman and Vasicek [22]). 

Our results are presented as follows. In Section [2] we present the two-factor spot model, and 
we compute analytical futures prices along with a discussion of pricing measures in Section 3. 
Section 4 explains in detail the estimation steps and the procedure applied to fit the model 
to data. The results of this estimation procedure applied to EEX data is presented and 
discussed in Section 5. We conclude in Section 6. 

Throughout we use the following notation. For a matrix D we denote by D* its transposed, 
and I is the identity matrix. For p 6 N we denote by e p the p-th unit vector. The matrix 
exponential e At is defined by its Taylor expansion e At = I + Y^=i with identity matrix 
/. We also denote by log + x = max(logx, 0) for i£l. 
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2 The spot price dynamics 



In most electricity markets, like the EEX, hourly spot prices for the delivery of 1 MW of 
electricity are quoted. As is usual in the literature on electricity spot price modeling, one 
assumes a continuous-time model and estimates it on the discretely observed daily average 
spot prices. We refer, for instance, to Lucia and Schwartz [2H] and Benth, Saltyte Benth and 
Koekebakker [3] for more details. 

We generalize the a-stable (C)ARMA model of Bernhardt, Kliippelberg and Meyer-Brandis [7] 
and Garcia, Kliippelberg and Miiller (23] by adding a non-stationary stochastic component 
in the trend of the spot dynamics. By modeling the trend as a combination of a stochastic 
process and a deterministic seasonality function rather than only a deterministic seasonality 
function, which seems common in most models, we are able to describe the low frequent 
variations of the spot dynamics quite precisely. As it turns out, this trend will explain a 
significant part of the futures price variations and lead to an accurate estimation of the risk 
premium in the EEX market. 

A two-factor spot price model for commodities, including a mean-reverting short-time dy- 
namics and a non-stationary long-term variations component was first suggested by Schwartz 
and Smith [35], and later applied to electricity markets by Lucia and Schwartz [25]. Their 
models were based on Brownian motion driven stochastic processes, more precisely, the sum 
of an Ornstein-Uhlenbeck (OU) process with a drifted Brownian motion. We significantly 
extend this model to include jump processes and higher-order memory effects in the dynam- 
ics. 

Let (fl, J 7 , {J-t}t>a, P) be a complete filtered probability space satisfying the usual conditions 
of completeness and right continuity. We assume the spot price dynamics 

S(t)=A(t)+Z(t)+Y(t), t>0, (2.1) 

where A is a deterministic trend/seasonality function and Z is a Levy process with zero 
mean. The process Z models the low-frequency non-stationary dynamics of the spot, and 
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can together with A be interpreted as the long-term factor for the spot price evolution. The 
process Y accounts for the stationary short-term variations. We will assume that Y and Z 
are independent processes. We follow Garcia et al. [23] and Bernhardt et al. [7] and suppose 
that Y is a stationary CARMA-process driven by an a-stable Levy process. 



2.1 The stable CARMA-process 

We introduce stationary CARMAQo, g)-Levy processes (see Brockwell [TT] ) and discuss its 
relevant properties. 

Definition 2.1 (CARMA(p, g)-Levy process). 

A CARMA(p, q)-Levy process {Y(t)}t>o (with < q < p) driven by a Levy-process L is 
defined as the solution of the state space equations 



with 



J p-2 



\ b P-l J 



Y{t) 
dX{t) 






v 1 ; 



b*X(t) 

AX(t)dt + e p dL(t), 



(2.2) 
(2.3) 



.4 








1 






1 



\ 





y ~ a p ~ a p~i ~ a p~2 



-ax J 



where ai, . . . , a p , 6q, . . . , fe p _i are possibly complex-valued coefficients such that b q = 1 and 
bj = for q < j < p. For p = 1 the matrix A is to be understood as A = —ai. 



The driving process L of Y will be a non-Gaussian a-stable Levy process {L(t)}t>o with 
characteristic function given by lnEe 12 ^'-* = tcfy^z) for z £ R, where, 



— 7° \z\ a 



1 — i/3(sign z) tan (™)) + i\iz for a ^ 1, 
7^1(1 + i/3^(sign z) log + ifiz for a = 1. 



(2.4) 
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The sign function is denned by sign z — — 1 for z < 0, sign z = 1 for z > and sign = 0, 
respectively. Further, a G (0,2) is the shape parameter, 7 > the scale, j3 G [—1,1] the 
skewness, and [i the location parameter. If 7 = 1 and \x = 0, then L is called standardized. 

The parameter a G (0, 2) determines the tail of the distribution function of L(t) for all t > 0. 
Moreover, only moments strictly less than a are finite, so that no second moment exists. 
This implies also that the autocorrelation function does not exist. For further properties 
on stable processes and Levy processes, we refer to the monographs by Samorodnitsky and 
Taqqu [33] and Sato 



The solution of the SDE (2.3) is a p-dimensional Ornstein-Uhlenbeck (OU) process given by 



X(t) = e A{t - s) X{s) + J e A{t - u) e p dL(u), < s < t, (2.5) 



where the stable integral is defined as in Ch. 3 of Samorodnitsky and Taqqu [33J. From (2.2) 
we find that Y is given by 

Y(t) = b*e Ait - s) X{s) + J b*e Ait - u) e p dL{u), 0<s<t. (2.6) 



Equations (2.2) and (2.3) constitute the state-space representation of the formal p-th. order 
SDE 

a(D)Y(t) = b(D)DL(t), t > 0, (2.7) 

where D denotes differentiation with respect to t, and 

a(z) := z p + a 1 z p - 1 + ■ ■ ■ + a p (2.8) 
b{z) := b + b 1 z + --- + b q z q (2.9) 



are the characteristic polynomials. Equation (2.7) is a natural continuous-time analogue of 



the linear difference equations, which define an ARMA process (cf. Brockwell and Davis [T3]). 

Throughout we assume that Y and X are stationary in the sense that all finite dimensional 
distributions are shift- invariant. Based on Proposition 2.2 of Garcia et al. [23] (which sum- 
marizes results by Brockwell and Lindner [15J) we make the following assumptions to ensure 
this: 
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Assumptions 2.2. Stationarity of CARMA-process. 



(i) The polynomials a(-) and b(-) defined in (2.8) and (2.9), resp., have no common zeros. 

(ii) E [log + |L(1)|] < oo. 

(Hi) All eigenvalues of A are distinct and have strictly negative real parts. 



Assumption (ii) and (iii) imply that X is a causal p-dimensional OU process, hence also Y 
is causal. 

Remark 2.3. Our model is a significant generalization of the two-factor dynamics by 
Schwartz and Smith [35] and Lucia and Schwartz [28]. Among various models, Lucia and 
Schwartz [28] suggested a two factor dynamics of the spot price evolution based on a short 
term Gaussian OU process and a long-term drifted Brownian motion. In our framework 
a Gaussian OU process would correspond to a Gaussian CARMA(1,0) process. It is clear 
that such a Gaussian model cannot capture the large fluctuations in the spot price, like for 
example spikes, and jump processes seem to be the natural extension. Based on the studies 
of Bernhardt et al. [7] and Garcia et al. [23], a-stable processes are particularly suitable 
for the short-term dynamics in the spot price evolution. Furthermore, empirical analysis of 
electricity spot price data from Singapore and Germany in [7J and [23] show strong statistical 
evidence for full CARMA processes to capture the dependency structure of the data. As 
for the long-term baseline trend, we shall see in Section [4] that a normal inverse Gaussian 
Levy process is preferable to a Gaussian process in a data study from the German electricity 
exchange EEX. □ 



2.2 Dimensionality of CARMA-processes 

A more standard model in electricity is to describe the spot by a sum of several OU processes, 
where some summands describe the spike behavior and others the baseline dynamics (see 
for example the model by Benth et al. [2]). A CARMA process is in a sense comparable to 
such models, as we now discuss. 
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By Assumption 2.2 iii) A has full rank, i.e. it is diagonalizable with e At = Ue Dt U~ 1 . Here, 
D is a diagonal matrix with the eigenvalues Ax, . . . , X p of A on the diagonal and U is the full 
rank matrix having the eigenvectors of A as columns. Since all eigenvalues have negative real 
parts, all components of e At are mean reverting. Each component of the vector e" 4 ^~^X(s) 



from (2.6) will therefore mean revert at its own speed, where the speed of mean reversion is 
a linear combination of the diagonal elements of e Dt . 



As we shall see in a simulation example of a CARMA(2,l)-process in Section 4.5, it captures 
the situation, where a first component has a slower rate of mean-reversion than the second 
(see Fig. pi). This is similar to a two- factor spot model, where the base and spike components 
of the spot price evolution are separated into two OU processes with different speeds of mean 
reversion. The advantage of working with a stable CARMA process, as we propose, is that 
it is possible to capture the distribution of the small and large jumps in one distribution. 
Since extreme spikes are rather infrequently observed, it is difficult to calibrate the spike 
component in a conventional two-factor model; this has been observed in Kliippelberg et 
al. [25]. With our CARMA- model, we avoid the difficult question of spike identification and 
filtering. 



3 The futures price dynamics 

In commodity markets, futures contracts are commonly traded on exchanges, including elec- 
tricity, gas, oil, and coal. In this section we derive the futures price dynamics based on the 



a-stable CARMA spot model (2.1). Appealing to general arbitrage theory (see e.g. Duffie 



[20] . Ch.l), we define the futures price f(t, r) at time t for a contract maturing at time r by 

f(t,r)=E Q [S(r)\jr t ], 0<t<r<oo, (3.1) 

where Q is a risk neutral probability measure. This definition is valid as long as S(t) G L 1 (Q). 
In the electricity market, the spot cannot be traded, and every Q ~ P will be a risk neutral 
probability (see Benth et al. [3]). For example, Q = P is a valid choice of a pricing measure. 
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In that case, the condition S(t) G L 1 (P) is equivalent to a tail parameter a of the stable 
process L being strictly larger than one, and a process Z with finite expectation. In real 
markets one expects a risk premium and hence it is natural to use a pricing measure Q ^ P. 



We will discuss possible choices of risk neutral probability measures Q in Section 3.1 



Based on our spot price model, we find the following explicit dynamics of the futures price 
for a given class of risk neutral probability measures: 



Theorem 3.1. Let S be the spot dynamics as in (2.1), and suppose that Q ~ P is such that 



L and Z are Levy processes under Q. Moreover, assume that the processes Z and L have 
finite first moments under Q. Then the futures price dynamics for < t < r is given by 

fit, r) = A(r) + Z(t) + b*e A ^X(t) + (r - t)E Q [Z{\)} + b* A' 1 (I - e A ^) e p E Q [L(l)] . 



Proof. Using f(t,r) = E Q [S(r) | T t ] = A(r) +E Q [Z(r) | T t ]+E Q [Y(r) | F t ). Since Z 

is a Levy process under Q, we find 

E Q [Z(r) | F t ] = Z(t) + E Q [Z(r) - Z(t) \ T t \ = Z(t) + (r - t)E Q \Z(\)\ . 

Now denote by M(u) = L(u) — Eq[L(1)]m for t < u < r, which has zero mean. Then, by 
partial integration, 



E, 



Q 



h * e A(r-u) ep 



U 



E Q [b*e p M(r) - b*e A(r -* ) e p M(0] - h* Ae A{T - u) e p E Q [M(u)\du = 



which implies E^ 



b*e A(T ~ u) e p dL( 



u 



E Q [L(1)) / b*e A ^e p du. 



Hence, the CARMA part of the spot dynamics converts to 



E Q [Y(r)\T t ]=E Q 



b*e A(T - 4) X(t) + / h*e A{T - u) e p dL(u)\F t 



b*e A(T -* ) X(t)+E Q 



b*e A{T - u) e p dL( 



u 



b*e A(T -* } X(t) + / b*e A(r - u) e p du E Q [L(l)] 



b*e A(T - 4) X(t) + WA- 1 (I - e A ^) e p E Q [L(l)] . 



Combining the terms yields the result. 



□ 
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In electricity markets the futures contracts deliver the underlying commodity over a period 
rather than at a fixed maturity time r. For instance, in the German electricity market 
contracts for delivery over a month, a quarter or a year, are traded. These futures are 
sometimes referred to as swaps, since during the delivery period a fixed (futures) price of 
energy is swapped against a floating (uncertain) spot price. The futures price is quoted as 
the price of 1 MWh of power and, therefore, it is settled against the average spot price over 
the delivery period. Hence, the futures price F(t,Ti,T 2 ) at time < t < T\ < T 2 for a 
contract with delivery period [Ti,T2] is defined as 

1 



F(t,T 1 ,T 2 )=E Q 



T 2 



S(r)dn 



(3.2) 



_T 2 —T\ J Tl 

where we have assumed that settlement of the contract takes place at the end of the delivery 
period, T 2 . 



Using Theorem 3.1 we derive by straightforward integration the swap price dynamics 



F(*,Ti,r 2 ) from (3.2) 



Corollary 3.2. Suppose all assumptions of Theorem 3.1 are satisfied. Then, 



F(t,T u T 2 ) 



T 2 — T\ J Tl 



Tz 



* /l-l 



A(r)dr + Z(t) 



b*A 



T 2 -T, 



(e A ^-e AT ^ e - At X(t) + T Q (t,T u T 2 ) 



where 



FQ(t,T 1 ,T 2 ) = Q(T 2 + T x ) - tj E Q [Z(1)} - |^ (e** - e^) e - At e p E Q [L(l)] 



+ b*A- 1 e p E Q [L(l)}. 



Proof. By the Fubini Theorem, we find 



F(*,T 1) T a ) = — ^— f 2 f(t,r)dr. 

J 2 — 1 1 JTi 



Applying Theorem 3.1 and integrating yield the desired result. 



□ 



The risk premium is defined as the difference between the futures price and the predicted 
spot, that is, in terms of electricity futures contracts, 



i? pr (t,T 1 ,T 2 ) = F(t,T 1 ,T 2 )-E 



T 2 



1 f T2 

— / S(r)dr\T t 

- 1 1 JTi 



(3.3) 
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From Cor. 



3.2 



we find that the theoretical risk premium for a given pricing measure Q is 



R(t,T h T 2 ) = TQ(t,Ti, T 2 ) — Tp(t, 7\, T 2 ) 

= Q(r 2 + ro-t)Eg[z(i)]-|^ 

+b*A- 1 e p (E Q [L(1)]-E[L(1)]) . 




)e~ At e p (Eg[L(l)]-E[L(l)]) 



(3.4) 



Here, we used the assumption that Z has zero mean under P. The first term gives a trend in 
"time to maturity" implied by the non-stationarity part Z in the spot price dynamics. "Time 
to maturity" is here interpreted as the time left to the middle of the delivery period. The 
two last terms are risk premia contributions from the CARMA short-term spot dynamics. 
They involve an explicit dependence on the speeds of mean-reversion of the autoregressive 
parts and the memory in the moving-average part. We will apply the risk premium in the 
empirical analysis of spot and futures data from the EEX. 

3.1 Equivalent measure transforms for Levy and a-stable pro- 



In this subsection we discuss a class of pricing measures that will be used for the specification 
of the futures price dynamics. 

We require from the pricing measure that Z and L preserve their Levy property and the 
independence. For this purpose, we consider probability measures Q = Ql x Q z , where Ql 
and Qz are measure changes for L and Z, respectively (leaving the other process unchanged). 
Provided Z has exponential moments, a standard choice of measure change is given by the 
Esscher transform (see Benth et al. [3], Section 4.1.1). Note that L, the a-stable process in 
the CARMA-dynamics, does not have exponential moments. 

We define the density process of the Radon-Nikodym derivative of Qz as 



cesses 



dQ z 



exp(e z Z(t)-<j> z (Oz)t), t>0, 



(3.5) 



dP F t 
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for a constant and <fiz being the log-moment generating function of Z(l) (sometimes 

called the cumulant function of Z). In order to make this density process well-defined, 
exponential integrability of the process Z up to the order of 9z must be assumed. Under 
this change of measure, the Levy measure of Z will be exponentially tilted by 9z, that is, 
if we denote the Levy measure of Z (under P) by u(dx), then its Levy measure under Qz 
becomes UQ z (dx) = exp(9 z )v(dx) (see Benth et al. [3], Section. 4.1.1-4.1.2 for details). 



To choose a risk neutral measure Ql is a more delicate task. We know from Sato [34J, 
Theorems 33.1 and 33.2, that equivalent measures Q exists for stable processes, however, it 
seems difficult to construct one which preserves the stable property. As an alternative, one 
may introduce the class of tempered stable processes (see e.g. Cont and Tankov [T7|, Chapter 
9), and apply standard Esscher transformation on these. 

A tempered stable process is a pure jump Levy process, where the stable-like behavior is 
preserved for the small jumps. However, the tails are tempered and, therefore, extreme 
spikes are less likely to be modeled with the tempered stable process. The Levy measure is 
given by 

c + e® LX c e® L \ x \ 

v TS {dx) = ^ 1+a l(o,oo)(g)dg+ i , 1+a l(-oo,o)(g)da;. (3.6) 

Here, 9l < and c_,c + G 1R+. A consequence of the tempering is that certain exponential 
moments exist. Tempering of a stable distribution results in a tempered stable distribution, 
and is analogous to taking an Esscher transform of the stable process using a negative 
parameter 6 L on the positive jumps, and a positive parameter — 6 L on the negative jumps. 

In particular, define g:li->las q(x) := c 9lX h 0jOO )(x) + e° L ^ lj-o^o)^) f° r some constant 
9l < 0. Suppose the stable distribution L has (under P) the characteristic triplet (jl, 0, Ul), 
where 

VL{dx) = l(o,oo)(aO dx + rrp^ l(-oo,o)(aO dx 
is the Levy measure of our stable process L. The parameters c + ,c_ can be matched to 



the parameters in (2.4), using Example 2.3.3 of Samorodnitsky and Taqqu [33]. Then the 
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tempered stable measure Ql, with characteristic triplet (7ts,0, Vts) is equivalent to the 
physical probability measure P (see Cont and Tankov [T7], Proposition 9.8), with drift 
parameter 7-^5 is given by 

{7l < a < 1 

(3.7) 
7l + f{\ x \<i} x(q(x) - l)v L {dx) 1 < a < 2 

and the Levy measure is given by VTs{dx) = q(x)vi,(dx). The special case of a Cauchy 
distribution (a = 1) is left out since one is not able to define the Levy-Kintchine formula using 
a truncation of the small jumps and the large jumps given by l| z |>i. For our applications it 
is of particular value to know the expectations of L(l) under P and Ql- 

Lemma 3.3. Let L be an a-stable Levy process under P with a G (1, 2). Find Ql by stable 



tempering for L < as in (3.6). Then the difference in mean of L(l) under Q L and P is 
given by 

E Qi [L(l)]-E[L(l)] = r(l-a)(-^r- 1 (c + -c_) , (3.8) 
where T is the gamma function. 



Proof. Using (3.7) and the Levy-Khintchine formula (e.g. Cont and Tankov [T7j, Prop. 3.13) 



for 1 < a < 2 we obtain 



E[L(1)) - E Ql [L(1)) = lL - its + [ x(v L - v TS ){d. 

J\\x\>l} 



X) 

{\x\>X} 

11 ~ Its + / x{l - q(x))v L {dx) 

J{\x\>l} 



c 



oo x JO 



= -T(l- a )(-e L ) a - 1 {c + -c) , (3.9) 

where we have used partial integration on the two integrals and rHospital's rule to obtain 
the last identity. This proves the result. □ 

Remark 3.4. By altering 8l one can match any relevant mean change in the risk premium 
Eq[L(1)] — E[L(1)], as long as this can be obtained by a negative choice of 6 L . This turns 
out to be appropriate for our applications. □ 
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4 Fitting the model to German electricity data 



Our data are daily spot and futures prices from July 1, 2002 to June 30, 2006 (available from 



http://eex.com). We fitted our model both to base load and peak load data, respectively. 
The futures contracts considered in this analysis are the Phelix-Base-Month-Futures and 
the Phelix-Peak-Month-Futures. Base load futures contracts are settled against the average 
of all hourly spot prices in the delivery period. Peak load futures contracts, on the other 
hand, are settled against the average of the hourly spot prices in peak periods of the delivery 
period. The peak period is counted as the hours between 8 a.m. and 8 p.m. every working 
day during the delivery period. The time series of daily spot prices used for our combined 
statistical analysis is taken to match the futures contracts: for the base load contracts we 
use the full time series consisting of daily observations including weekends (i.e., we have 7 
observations per week), while in the case of peak load contracts the weekends are excluded 
(i.e., we have 5 observations per week). 

Figures [T] and [2] show the spot and futures prices for both base load and peak load. From 
these plots we can see similar patterns of the base and peak load data, however, peak load 
data are more extreme. Note that all plots cover the same time period; however, for the 
base spot data we have 1461 observations, whereas for the peak spot data we have only 1045 
observations in the same period, due to the missing weekends. 

The estimation procedure for our model consists of several steps, which are explained in the 
following. 



4.1 Seasonality function A 

The estimation of the deterministic trend component A is a delicate question. A mis- 
specification of the trend has a significant effect on the subsequent analysis, in particular, on 
the risk premium. Motivated by the seasonality functions used in Bernhardt et al. [7] and 
Garcia et al. [23] , we take the seasonality function of the peak load contracts as a combination 
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Figure 1: Daily spot prices from July 1, 2002 to June 30, 2006, base load (top) and peak load 
(bottom). 

of a linear trend and some periodic function 



Note that we choose a slightly simpler seasonality function than Bernhardt et al. [7j and 
Garcia et al. [23J, only taking the mean level, a linear trend and a yearly periodicity (modeling 
the weather difference between summer and winter) into account. Weekly periodicity in peak 
load contracts is not that pronounced, since weekends are not considered in peak load data. 
Since no new trading information is entering during the weekend (trading takes place during 
weekdays), we will adjust the periodicity to 261 and consider the peak load contracts as a 
continuous process on all non-weekend days. 

In the base load prices a clear weekly seasonality is visible. Weekend prices are in general 



A p (t) 
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future prices in Euro (base) 
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Figure 2: Daily futures prices from July 1, 2002 to June 30, 2006, base load (top) and peak 
load (bottom). 



lower than during the rest of the week, and over the week one observes the pattern that 
Monday and Friday have prices lower than in the middle of the week. Therefore we include 
a weekly term in the base load seasonality function: 

. . /2ttA /2ttA /2ttA /2ttA , , 

A 6 (t) = ci + c 2 t + c 3 cos I— J +c 4 sm I — I + c 5 cos I — I +c 6 sm I — I .(4.2) 

Since time t is running through the weekends, a yearly periodicity of 365 is chosen. 

In the following, we will analyse both data sets, base load data and peak load data (spot 
and futures, respectively). For simplicity we will suppress the indices p for peak load and b 
for base load. 

This seasonality functions can be estimated using a robust least-squares estimate on the 



data; see Table 4.1 for the resulting estimates. The first plot of Figure M shows the 1461 
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Cl 


c 2 


c 3 


c 4 


c 5 


c 6 


Base 


19.4859 


0.0217 


-2.8588 


0.6386 


-6.7867 


2.8051 


Peak 


30.7642 


0.0349 


-2.5748 


1.5762 







Table 4.1: Estimated parameters of the seasonality function A(-). 



observations of the base data set together with the estimated trend and seasonality curve, 
the second plot zooms into this plot and shows these curves for the first 200 observations. 
From this second plot we can clearly see that the daily seasonality is caputered by A quite 
accurately. For the peak data over the same period (consisting of 1045 observations) we 
get similar plots (hence, we omit them here), except of the fact that A does not show a 
weekly seasonality in this case. The overall growth rate was very small during the data 
period, justifying a linear term as a first order approximation. Note that we have introduced 
the non-stationary stochastic process Z to absorb all stochastic small term effects in the 
seasonality. This term will play a prominent role for the futures prices later. 

Subtracting the estimated seasonality function from the spot data leaves us with the reduced 
model Z(-) + Y(-), where we have neglected in our notation the fact that we have subtracted 
only an estimator of A(-). 

Next we want to estimate both components Z and Y invoking the deseasonalised spot price 
data and the futures prices. We will exploit the fact that the futures prices far from delivery 
will have a dynamics approximately given by the non-stationary trend component Z. Only 
relatively close to delivery, large fluctuations in the spot price dynamics are reflected in the 
futures prices. Since it is not clear how far away from delivery we need to be before the 
approximation of futures prices by Z works well, we will invoke an optimization routine to 
find the optimal distance. For this purpose, we introduce the notation u := |(Ti + T 2 ) — t, 
which will be referred to as "time to maturity". 

Denote by u* the optimal time to maturity (we will define what we understand by "optimal" 
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Figure 3: Base spot prices and estimated seasonality function. Top: whole period (14-61 
observations) . Bottom: first 200 observations. 



below), where futures contracts with time to maturity u > u* have a dynamics approximately 
behaving like the non-stationary term. How big to choose it* is not possible to determine a 
priori, since we must analyse the error in an asymptotic consideration of the futures prices 



see (4.3) and (4.4) below). This error is highly dependent on the parameters in the spot 



price model, which we do not yet know. In the end, u* should be chosen so that the error in 



the risk premium estimation is minimal; cf. Section 4.6 



The estimation will, therefore, be repeated using different values of u* (this parameter will 
sometimes be called threshold m. the following); i.e., we choose a subset U* := \u ^ n , «max] — 
[v/2,Mf], where v is the average delivery period and Mf is the maximal time to maturity 



observed in the futures data set, and perform the steps of the Sections 4.2-4.6 below repeat- 



edly for u* & U* . For each value u* G U* the error in the risk premium is calculated, and u* 
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is the value which minimizes this error among all u* G U*. This optimal threshold u* is then 
considered as final choice of u* for the calculation of all estimates including the processes Z 
and Y and for the CARMA parameters. 

One should keep in mind that for too large u* there is only few data available with u > u*, 
which yields unreliable estimates. Since we count the time to maturity as number of trading 
days until the mid of the delivery period (which has length v), time to maturity is always at 
least v/2; hence we do not consider any u* smaller than v/2. Overall, we decided to choose 
u*: = \v/2] (which is 16 for the base load and 11 for the peak load data) and «max = Mf/2 
(note that Mf is 200 days for the base load and 144 days for the peak load contracts). As 
we will see later, the optimal u* in our data examples is quite small, so that this choice of 
li^ax i s completely satisfying. 

Next we want to explain in detail, how we separate Z and Y for a given fixed time to 
maturity. Consequently, we perform the model estimation for all ^(7\ + T 2 ) < u* < 200(146) 
and take all futures prices for the estimation procedure, whose time to maturity u > u* . 

We explain each step in the estimation procedure in detail: 



4.2 Filtering the realization of the non-stationary stochastic pro- 
cess Z 



Recall the futures price F(t, T\, T%) in Corollary 3.2 Since we assume that the high-frequency 
CARMA term Y is stationary, it holds for fixed length of delivery T 2 — T\ that 

lim (e** - e**) e-*X(f) = (4.3) 

lim ^^(e^-e^)e-^e P E Q [L(l)] = 0. (4.4) 

Hence, in the long end of the futures market, the contribution from Y to the futures prices 
may be considered as negligible. In particular, from the futures price dynamics (Corol- 
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lary 3.2) we find for [Ti,T2] far into the future (that is, t much smaller than T\) that 

F{t,T u T 2 ) := F(t,n,T 2 ) - -±— f\{r)dr 

^Z(t) + b*A- 1 e p E Q [L(l)]+ Q(T 1 + T 2 )-^ E Q [Z(1)). (4.5) 

Recalling the notation u := |(Ti + T 2 ) — t coined "time-to-maturity", we slightly abuse the 
notation and introduce F(t,u) := F(t,Ti,T2). 

For u > u*, we approximate 

H P (u) :=E[F{t,u)\ 

« E[Z(i)] +bM- 1 e p E Q [L(l)] + «E Q [Z(1)] 
= bM- 1 e p E [I(l)] + t 1 E Q [Z(l)] 

=: C + uE Q [Z(l)], (4.6) 

where we have used the zero-mean assumption of Z under P. This approximative identity 
can now be used for a robust linear regression on the time to maturity u, in order to estimate 
the real numbers C and Eq[Z(1)]. Knowing these two parameters enables us to filter out 



the realization of the process Z. According to Equation (4.5) we obtain 
Z(t) = Z (\{ Tl + T 2 ) -u)= * E [F(t, T U T 2 ) - C - u%[Z{l)\ 

(4.7) 

where U(t,u*) := {{u,T u T 2 ) G M 3 | u > u* and 3F(t,T 1 ,T 2 ) : \{T l +T 2 ) - t = u). 

Remark 4.1. Note that after estimating the CARMA parameters, we can also find an 
estimate for E Q [L(1)] simply by taking E Q [L(1)] = C(b*A -1 ep) -1 . □ 

Remark 4.2. We recall that the futures market at EEX is not open for trade during the 
weekend. Therefore, using our estimation procedure, we do not get any observations of Z 
during weekends. We will assume that Z is constant and equal to the Friday value over the 
weekend, when filtering the non-stationary part of the spot in the base load model. One 
may argue that this strategy could lead to large observed jump of Z on Monday morning, 
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when all information accumulated over the weekend is subsumed at once. We will return to 
this question in Section [5\Tj □ 



4.3 Estimation of the CARMA parameters 



Recall our spot-model (2.1 ) 



S(t)=A(t) + Y(t) + Z(t), t>0. 



After A(-) and Z(-) have been estimated in Subsections 4.1 and 4.2, respectively, a realization 
of the CARMA-process Y can be found by subtracting both from the spot price. Figure 
[3] shows the estimated processes Z (dotted red) and Y + Z (black) for both the base and 
the peak load data, for the full period July 1, 2002 to June 30, 2006 and for u* = 16 
exemplarily, each. Obviously, the process Z captures the medium-range fluctuations and Y 
the short-range fluctuations of the detrended and deseasonalized process Y + Z. 

Again we keep in mind that the process Y is the result of some estimation procedure. There 
exists a number of papers devoted to the estimation of the CARMA parameters in L 2 (see 
for instance Brockwell et al. [H], Tsai and Chan [36]). Methods can be based either directly 
on the continuous-time process or on a discretised version. The latter relates the continuous- 
time dynamics to a discrete time ARMA process. The advantage of this method is obvious, 
since standard packages for the estimation of ARMA processes may be used in order to 
estimate the parameters of the corresponding CARMA process. Some care, however, is 
needed since this approach does not work in all cases. Brockwell and collaborators devote 
several papers to the embedding of ARMA processes in a CARMA process; cf . pUl EE] . Not 
every ARMAQo, q) process is embeddable in a CARMA (p, q) process. 



From Assumptions 2.2, it follows by Proposition 2 of Brockwell et al. [H] (cf. also Garcia 
et al. [23], Prop. 2.5) that every CARMA(p, q) process Y observed at discrete times can be 
represented as an autoregressive process of order p with a more complex structure than a 
moving average process for the noise. For such a discretely observed CARMA process Y on 
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Figure 4: Estimated processes Z (red) and Y + Z (black) in the period July 1, 2002, to June 
30, 2006, for both the base load data (top) and the peak load data (bottom). 

a grid with grid size h, denoting the sequence of observations by {y n } n ^; i.e. y n = Y(nti), 
Prop. 3 of Brockwell et al. [H] gives 

p 

- e Xi B)y n = e n . (4.8) 

i=l 

Here, B is the usual backshift operator and {£ n }neN is the noise process, which has repre- 
sentation 

P [-nil 

s n = J2 Ki II^ 1 - / e^ {nh - u) dL{u) . (4.9) 

i=l j^i J(;n-l)h 

The constants are given by «j := 6(Ai)/a'(Aj) and Ai,...,A p are the eigenvalues of A. 
The process {e n }neN is p-dependent. When L has finite variance, e n has a moving average 
representation; cf. Brockwell et al. [Mj, Proposition 3.2.1. However, for the case of infinite 
variance this is no longer true and causes problems. 
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Davis [T8] . Davis and Resnick [19] and Mikosch et al. [29] have proved that ordinary L 2 - 
based estimation methods for ARMA parameters may be used for a-stable ARMA processes, 
although they have no finite second moments. Moreover, Davis and Resnick [T9] showed 
that the empirical autocorrelation function of an a-stable ARMA process yields a consistent 
estimator of the linear filter of the model, although the autocorrelation function of the 
process does not exist. Hence the parameters ax,...,a p can be consistently estimated by 
L 2 -methods. In [19J it has also been shown that the rate of convergence is faster than in the 
L 2 -case. 

Already in the analysis performed in Garcia et al. [23] the CARMA(2,1) process has been 
found to be optimal. Although our model is slightly different, it turns out that this CARMA 
dynamics is still preferrable for Y based on the AIC model selection criterion. Hence, in the 
following example we spell out the above equations for the case of a stable CARMA(2,1) 
model. 

Example 4.3. [The CARMA(2,1) process] 

By applying (4.8) and ( 4.9[ ) for the case of a CARMA(2,1) process, we find the discrete-time 



representation for a gridsize h > 0, 

y n = (e Xlh + e X2h ) y n ^ - e^ +x ^ h y n _ 2 + e n , 

where e n is given by 

rnh 

e n = ( Kl e x ^ nh - u) + K 2 e X2{nh ~ u) )dL{u) 

J(n-l)h 
p(n-l)h 

+ / {K ie X2h e Xl{nh - u) + K 2 e Xlh e X2{nh - u) )dL(u). 

J(n-2)h 

The two integrals in the noise are independent. It is, however, not possible to recover the 
noise by simple multiplication and subtraction as in the ARMA case. The actual relation of 
two successive noise terms e n and e n+ i is based on the continuous realization of {L(t)} t >o in 
the relevant intervals, which is unobservable. □ 

For the mapping of the estimated ARMA parameters to the corresponding CARMA param- 
eters we observe that equation (4.8) is a complex way to express that {e~ Xih } i= i^ ^ p are the 

26 



roots of the autoregressive polynomial 4>(z) = 1 — 4>iz — • • • — <p p z p of the ARMA process. 
We proceed, therefore, as follows for identifying the CARMA parameters from the estimated 
ARMA process: 

• Estimate the coefficients 0i, . . . , P of the ARMA process 

• Determine the distinct roots £j for % = 1, . . . ,p of the characteristic polynomial. 

• Set Xi = — log(£j)//i, where we recall that h denotes the grid size. 

Because of the simple structure of the autoregressive matrix A of the CARMA process we 
can calculate the characteristic polynomial P of the matrix A as 

P(A) = (A p + a^- 1 + • • • + op) . 

Since the Aj are the eigenvalues of A, we know that -P(Aj) = for i — 1, . . . , p. Hence, given 
the eigenvalues Aj of the matrix A we recover the coefficients a±, . . . , a p by solving a system 
of p linear equations. 

We estimate the moving average parameters based on the autocorrelation function. For 
its estimation we apply a least absolute deviation algorithm based on the empirical and 
theoretical autocorrelation functions of the CARMA process. The theoretical autocorrelation 
function of y takes the form 

7 tf (s) = b*e A|s| £b, s > 0, 

where the matrix S is given by 

POO 

S = / e Au e p e* p e A * u du = -A^e^; . 
Jo 

In this representation, A -1 is the inverse of the operator A:X4 AX + XA* and can be 
represented as vec -1 o ((A ® I p ) + (I p <g> A))' 1 o vec (see Pigorsch and Stelzer [3 1J ) . Using the 
above procedure for the estimation of the moving average parameter b is based on second 
order structure and, therefore, not straightforward to use for stable processes. In practice 
this procedure works and we can use it to estimate the moving average parameter b. 
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4.4 Estimation of the stable parameters 



After estimating the autoregressive parameters, the noise {e n } ne N can be recovered. Recall 
that, by p-dependence, the noise terms of lags m > p are independent. Motivated by results 
for discrete-time stable ARMA processes, Garcia et al. [23] have applied estimation methods 
for independent noise variables. They have also shown in a simulation study that one gets 
quite reliable estimates by treating R := {e„} ne N as independent sequence. 

By a simple computation we can relate the estimated parameters of the series Ri back to an 
estimate of the a-stable process L. We show this for the CARMA(2,1) model in the next 
example: 



Example 4.4. [Continuation of Example 4.3 cf. [23] 



Using Samorodnitsky and Taqqu [33J, Property 3.2.2, a relation in distribution between 
the a-stable process L and the noise process {£ n }neN of the ARMA(2,1) model sampled on 
a grid with grid size h can be established (for clarification, we will now sometimes write 
Ili 0l, A*l) for the parameters of the a-stable process L, which were so far denoted just 
by (a, 7, (3, /i)). In particular, e n has an a-stable distribution with parameters (a e , j E , f3 £ , // e ) 
given by 

a £ = (Xl = ol 

le=U |«i e Xl{h ~ u) + k 2 e X2{h ~ u) | Q + | «4 e X2h e Xl{h - u) + k 2 e Xlh e X2{h ~ u) |° dii) ^-/l 



l^e — A*L — A* fo r ^7^1 

Note that, for a and p being real numbers, a' p ' := |a| p sign(a) denotes the signed power 
(Samorodnitsky and Taqqu [33], eq. (2.7.1)). Moreover, we can easily see that (3 e = 0l, if 
both Ki and k 2 are positive. 



h 



28 



4.5 Recovering the states 



In order to calculate the theoretical futures prices derived in Corollary 3.2 it is necessary to 
recover the states X of the CARMA-process. Brockwell et al. [H] describe a rather ad-hoc 
method to do this by using an Euler approximation. 



In the linear state space model (2.2), the Kalman filter is the best linear predictor provided 
the driving noise is in L 2 . Since a-stable Levy processes for a G (0, 2) do not have finite 
second moments, the Kalman filter will perform unsatisfactorily. One possibility to resolve 
this is to apply a particle filter, which does not require a finite second moment of the noise 
process. However, the particle filter requires a density function instead, which poses a new 
problem for a-stable processes. Integral approximations of a-stable densities exist, but they 
are time consuming to calculate and simple expressions do not exist. One can use a particle 
filter by simulating from the a-stable distribution, but this is also very time consuming. A 
large number of paths need to be simulated in order to get a reasonable estimation (even when 
using appropriate variance reducing methods like importance sampling). As an attractive 
alternative, we introduce a simple L 1 -filter applicable to CARMA processes with finite mean. 



Recall from (2.2)-(2.6) that we can work with the following state-space representation of the 



CARMA process 



Vn = b*x n , (4.10) 



"nh 

x„ = e A/l x n -i + z n with z n = / e A{nh - u) e p dL(u) (4.11) 

J(n-l)h 

Here, y n and x n are discrete observations of Y and X, respectively, on a grid with grid size 
h. 

Notice that given y n and x n _x the value of b*z n is determined and given by 

b*z n \y n , x n _i =y n - b*e A/l x n _i. (4.12) 
This will come to use in a moment when deriving the filter. 
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First, we make an "Euler" approximation of the stochastic integral defining z n by 



z„ 



h 



nh 



(n-l)h 



jA(nh—u) 



e p duAL{n, h) = -A" 1 (I - e Ah ) 



AL(n,h) 
h 



where AL(n,h) = L(nh) — L((n — l)h). Note that a traditional Euler approximation (see 
Kloeden and Platen (23]) would use the left end-point value of the integrand in the approxi- 
mation, whereas here we use the average value of the integrand over the integration interval. 
We find 



E[z n | y n , x n _i] « -E[AL(n, h)/h \ y n , x n _ x ] A 1 (/ - e Ah ) e p . 



(4.13) 



Multiplying (4.13) with b* and combining it with (4.12) gives 



E[AL(n, h)/h | y n ,x n _i] 



y n — b*e Afe x. 



n-l 



-b*A- 1 (I - e Ah ) e p ' 



By plugging (4.14) into (4.13) we find 

E[z n | y n , x n _x] w -A" 1 (/ - e Ah ) e p 



y n - b*e Afe x n _x 
-bM- 1 (J-e A/l )e p ' 



(4.14) 



(4.15) 



We can use this as an L 1 -filter for z n . Applying (4.15), we can filter the states X of the 



CARMA-process. Using the state equation (4.11) we find 



E[x n | y n , Xn_i] w e j4/l x„_ 1 + E[z n | y n , x„,_!]. 



(4.16) 



We tested the filter on simulated data from a CARMA(2,1) process with the same param- 



eters as we find from our model for the base spot prices (see Table 5.3). The path of the 
CARMA(2,1) process was simulated based on an Euler scheme on a grid size of 0.01 for 
< t < 1461, and the a-stable Levy process was simulated using the algorithm suggested 
by Chambers, Mallows and Stuck [16J. The estimation of the states is done on a grid with 
grid size h — 1. Figure [5] shows the estimated states (red curve) for both state components 
together with the simulated states (black curve). It is clearly visible that the L 1 -filter gives 
a good approximation of the true states X driving the a-stable CARMA process Y. 
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simulated and filtered states (first component) 
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Figure 5: Estimated states (red) and true states (black) of a simulated CARMA(2, 1) process 
using the ^-filter. 

4.6 Risk premium comparison 

In order to find the optimal threshold u* for filtering out the non-stationary process Z from 
the futures data, we compare the empirically observed risk premium with its theoretical 
counterpart. 



Recall the risk premium R pr in (3.4) implied by the futures price dynamics in Corollary 3.2 
By using v = T<i — T\ and recalling the notation u = \{Ti + T%) — t, we can rewrite R pr for 
u > 1(^2 — T\) and fixed v (being one month in our studies) to 

R VT (u\u,v) = ~\h*A- 2 (e~* Av - e"^ e Au e p (E Q [L(1)] - E[L(1)]) 

+ b*A- 1 e p (Eq[L(1)]-E[L(1)])+uE q [Z(1)]. (4.17) 



Note that we can estimate all parameters in (4.17) only depending on a chosen threshold 
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u* . Hence, the risk premium in Equation (4.17) also depends on u* . In order to find an 



optimal u* we compare the risk premium (4.17), which has been estimated on our model 



assumptions, with the mean empirical risk premium based on the futures prices given by 
R pr (u*,u,v) := ^b*A- 2 (e^ Av - e~* Av> ) e A "e p E[L(l)] - bM- 1 e p E[L(l)] 



— : > \F(t,T u T 2 ) - f 2 K{r)dr (4.18) 

L 1 2 — 1 1 JTi 



b " A ~' (e AT2 -e ATl )e- At X(t)-Z(t) 



T 2 -T x 
Here, 

U(u, v) := |t, Tx, T 2 eR;^(T 2 + T 1 )-t = u,T 2 -T 1 = v and F(t, T lt T 2 ) exists } . 

The dependence on the threshold u* is only implicit. The estimated sample paths of Z 
and Y depend on u*, therefore, also the CARMA parameters A, b, the stable parameters 
(a,/3,7,/x) and the estimated sample paths of the states X also depend on u*. In order 
to compute R pr and R pr all these estimated parameters are used. Consequently, by using 
different thresholds we will get different estimates and different risk premia. We want to 
choose an optimal threshold u*, such that the mean empirical risk premium R rp is as close 
as possible to the model based risk-premium R pr . We invoke a least squares method, i.e. we 
minimize (for fixed v) the mean square error between the two functions (R pr (u*, u, v )) u> i v 
and (Rpr(u*,u,v)) u> i v with respect to all chosen thresholds u* G U* for the estimation 



procedure, cf. Section 4.1 



M f 

X It-./* \ t~* / =k \|2 



u* = argmin u , gC/ , >^ \R pr (u*, u, v) — R pr {u* : , u, v) 

u=v/2 

Here the dependence of the error function 

f(u*,v) := \R P r(u*,u,v) - R pr (u*,u,v)\ 2 (u*EU*) (4.19) 

u=v/2 

on u* is only implicit. In our data v corresponds to the average number of days per month 
(i.e. v = 1461/48 = 30.44 for the base data and v = 1045/48 = 21.77 for the peak data), and 
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the number Mf is the longest time to maturity, which we recall from Section 4.1 being 200 



for the base load contracts and 144 for the peak load. In order to calculate this minimum 
we calculate the values of f(u*,v) for all u* G [u/2, M//2] fl N. Figure [6] shows the risk 
premium error function f(u*,v) for base load (left) and peak load (right). In both cases, 
the minimum is attained at u* = 16. So our estimation procedure considers only base load 
forward contracts with delivery at least (about) two weeks away, and peak load contracts 
with delivery at least (about) three weeks away. 



error function f(u*,30.44) (base) 



20 



40 



60 



80 



100 



error function f(u*,21 .77) (peak) 



10 20 30 40 50 60 70 



Figure 6: The risk premium error function for base data (left) and peak data (right), for 
v/2 < u* < Mf/2. In both cases, f has its minimum at u* — 16. 



Below we present a summary of the estimation algorithm again. 



The algorithm 



Estimate A(-) as in (4.2) for the base load model and as in (4.1) for peak load model, and 
subtract from S(-). 



For each threshold u* G U*\ 



Approximate fip(u) = C + mEq[Z(1)] for u > u* and estimate C, Eq[Z(1)] by linear 



regression (4.6) 
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filter Z by flO; 



model Y = S — A — Z as CARMA(2,1) process, estimate the coefficients a\, 0,2, bo (recall 
that bi = 1) and estimate the parameters jl, Pl, Hl) of L; 



estimate Eg[L(l)] using (4.7) 



filter states of X = (X h X 2 )* using (4.16) 



calculate R pr (u*,u,v) as in (4.17) invoking the estimated parameters and states from 
the former steps; 



calculate R pr (u*,u,v) as in (4.18) invoking the estimated parameters and states from 



the former steps and the futures data. 



Now define the mean square error of the estimated R pr (u* ,u,v) and R pr (u*,u,v) based on 
all different thresholds u* G U* . The optimal threshold is found to be u* . 



5 Estimation results 



We now report the other results from the estimation procedure, when using the the optimal 
threshold u* = 16 both for base and peak load data. In this section we discuss the estimated 
values and their implications. 



5.1 Distributional properties of the filtered sample path of Z 



For the filtered Z which was found using (4.7) we can derive certain properties. Both for the 



base and the peak data, the realization of Z shows uncorrelated increments, see Figure [7| 

Figure [8] shows QQ-plots for the increments of Z versus a corresponding normal distribution, 
both for base data (left) and peak data (right). Note that the empirical variances of the 
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acf for increments of Z (base) acf for increments of Z (peak) 



Figure 7: Empirical autocorrelation functions for the increments of Z , for base data (left) 
and peak data (right). 



increments are 0.35 and 2.78 for the base and peak data, respectively. From these plots we 
conclude that for both data sets the increments of Z have heavier tails than the Gaussian 
distribution, and that this feature is even more pronounced for the peak data. 



QQ-plot for increments of Z (base) QQ-plot for increments of Z (peak) 
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quantiles of N(0,0.35) quantiles of N(0,2.7S) 



Figure 8: QQ-plots for the increments of Z against suitable normal distributions, for base 
data (left) and peak data (right). The reference distributions are N (0,0.35) and N (0,2.78) 
for base and peak data, respectively. 



Kernel density estimates suggest that the increments of Z can be described quite well using a 
normal inverse Gaussian (NIG) distribution (we refer to Barndorff-Nielsen [1] for a thorough 
discussion on the NIG distribution and its properties). The red curves in Figure M show log- 
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increments of Z: log-density (base) increments of Z: log-density (peak) 




-2-1012 -50 5 



Figure 9: Log-densities for the incremets of Z (red) as well as for NIG distributions (black 
solid) and normal distributions (black dashed) that have been fitted to these increments, for 
base data (left) and peak data (right). 





a z 


Pz 


5z 




Base load 


0.6451 


0.0998 


0.2206 


-0.0346 


Peak load 


0.2371 


-0.0083 


0.6582 


0.0230 



Table 5.2: Estimated parameters of the NIG distribution Q for the increments of Z. Since we 
assume that K[Z] = 0, the parameters have been estimated conditionally on this assumption. 



density estimates for the increments of Z for the base data (left) and the peak data (right), 
respectively. For comparison, we also plot the log-density curves of NIG distributions (black 
solid curves) and normal distributions (black dashed curves) that have been fitted to the 
increments of Z via maximum likelihood (the parameters for the NIG distributions can be 



found in Table 5.2). Clearly, the NIG distribution gives a much better fit than the normal 
distribution. Hence, we identify the non-stationary process Z with a normal inverse Gaussian 
Levy process. 

Remark 5.1. Recall that futures contracts are only traded on weekdays and, therefore, 
no variability in Z during weekends is observed. For base load contracts we have thus 
assumed that Z is constant during weekends for filtering purposes, but only considered 
weekdays data for analysing the distributional properties of Z. As we already mentioned in 



Section [4721 this strategy could lead to larger up- or downward movements of Z on Mondays, 
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when all information accumulated over the weekend is subsumed at once. However, we do 
not find such a behaviour in the estimated increments of Z . To see this, we calculated 
the variance of the increments of Z in the base data set for different weekdays separately. 
We find a variance of 0.3723 for the increments occurring from Fridays to Mondays, and a 
variance of 0.3526 for the remaining increments, i.e. Mon/Tue, Tue/Wed, Wed/Thu, Thu/Fri. 
For comparison, the corresponding variance from Wednesdays to Thursdays only is 0.3202, 
whereas for the remaining increments, i.e. Mon/Tue, Tue/Wed, Thu/Fri, Fri/Mon, the value 
is 0.3649. Furthermore, the estimated overall variance of the increments of Z is 0.35. Taking 
the estimation error for all these variances into account, we do not observe a significantly 
higher variance for the increments from Friday to Monday in the base data. This supports 
the idea that no relevant information enters the futures market over the weekends. Maybe 
this even can be expected, since the futures deal with a product to be delivered quite far 
in the future; hence, only really influential news with a long-range impact should affect the 
market. 

□ 



5.2 Estimation of the CARMA parameters 

The autocorrelation and partial autocorrelation function of the data suggest that there are 
two significant autoregressive lags, but also a relevant moving average component. Also 
the AIC and BIC criterion confirm that a CARMA(2,1) model leads to the best fit. The 



estimated parameters of a CARMA(2,1) model are given in Table 5.3 



For the estimated parameters (01,02) in the autoregression matrix A the eigenvalues of A 
are real and strictly negative, being Ai = —0.0641, A2 = —1.4213 for the base load and 
Ai = —0.1014, A2 = —2.2319 for the peak load. Our parameters satisfy Assumptions 2.2 
Hence, the estimated model is stationary. 

The estimates of oll (1.6524 for base and 1.3206 for peak) confirm that extreme spikes are 
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CARMA parameters 


Stable parameters 


a x 


a 2 


b 


a L 


Pl 


IL 


E[L(1)] 


Base load 


1.4854 


0.0911 


0.2861 


1.6524 


0.3911 


6.4072 


0.0566 


Peak load 


2.3335 


0.2263 


0.6127 


1.3206 


0.0652 


6.5199 


-0.0448 



Table 5.3: Estimates of the CARMA parameters and of the parameters of the stable process 
L. 







C 




0z 




Base load 


-0.0243 


1.6587 


-0.5282 


-0.1093 


— 0.0021 


Peak load 


-0.0382 


3.5678 


-1.3178 


-0.0168 


— 0.0552 



Table 5.4: Estimates of parameters determining the risk. 



more likely in the peak load data. As we can conclude from the positive signs of the skewness 
parameter /3l, positive spikes are more likely to happen than negative spikes for both data 
sets. Note, however, that a direct comparison of the values of /3l for the base and the 
peak load data is misleading, due to the significantly different parameters a^. Indeed, if 



we calculate the empirical skewness for the estimated e n (cf. Section 4.4) directly, we get a 
value of 0.28 for the base load data and 1.59 for the peak load data (cf. the comment on 



f3 £ = 0i in Example 4.4 for the base load data (/%, k 2 ) = (0.1636, 0.8364), and for the peak 
load data («i,k 2 ) = (0.2400,0.7600)). 



5.3 Market price of risk and risk premium 

We next present the results on the risk premium and the parameters for the market price of 
risk, based on our statistical analysis of base and peak load contracts with threshold u* = 16 



days in both cases. Estimates of the relevant parameters are presented in Table 5.4 



Recall that Eq[Z(1)] and C are found from regression (4.6); using the estimates of the 
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CARMA model from Table 5.3, we derive an estimate of Eg[L(l)]. Having both Eq[L(1)] and 
Eq[Z(1)] estimated, we can compute the parameters in the respective measure transforms 
of the NIG Levy process Z and the stable process L. For an NIG Levy process we use 
the fact that an Esscher transformed NIG(az, 0z,8z, \iz) random variable Z is again NIG 
distributed with parameters (atz, fiz + &zi $z, Hz) (see e.g. Benth et al. [3], p. 99). Using the 
mean of an NIG distributed random variable it holds that 



E Q [Z(l)] = fi z + 



V<Xz-Wz + 0z) 2 ' 

where 9z is the market price of risk for Z. Since estimates for the parameters az, 0z, &z 



and pLz are known from Table 5.2 we can use the estimate for Eq[Z(1)] together with the 
above equality to obtain an estimate of 6z, which results here in 6z = —0.1093 for the base 
load and in 6z = —0.0168 for the peak load data. Since 6z is estimated negative, more 
emphasis is given to the negative jumps and less emphasis to the positive jumps of Z in the 
risk neutral world Q. We see from the estimate on the risk-neutral expectation of Z that 
the contribution from the non-stationarity factor of the spot on the overall risk premium is 
negative. This is natural from the point of view of the hedging needs of producers. The 
non-stationary factor induces a long-term risk, which is the risk producers want to hedge 
using futures contracts. 

Lucia and Schwartz [28] also find a negative market price of risk associated to the non- 
stationary term in their two-factor models, when analysing data from the NordPool market. 
We recall that they propose a two-factor model, where the non-stationary term is a drifted 
Brownian motion. The negative market price of risk appears as a negative risk-neutral drift, 
which corresponds to a contribution to the risk premium similar to our model. We refer 
to Benth and Sgarra [5] for a theoretical analysis of the Esscher transform in factor models 
applied to power markets. 

Using the relations 7° = c + + c_ and (3 = (e+ — c_)/(c+ + c_) in the stable parameters as 
calculated in Example 2.3.3 of Samorodnitsky and Taqqu [33] and plugging in the estimated 
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parameters from Table 5.3, we find 



c+ = hi + /3 L ) 7 r = 14.9715 and c_ = hi - L )ff = 6.5532 



for the base load data, and 



c+ = i(l + /3 L ) 7 r = 6.3342 and c_ = 1(1 - fofr^ = 5.5587 



for the peak load data. Then by using (3.9) we can derive an estimate for 9l by 



0/ 



Eq[L(l)]-E[L(l)] 
T(l - al)(c + - c_) 



which leads to 8l = —0.0021 for the base lead data and to 6i = —0.0552 for the peak load 
data. The market price of risk for the CARMA-factor noise L is also negative, however, unlike 
the non-stationary factor a negative sign does not necessarily lead to a negative contribution 
to the risk premium. As we already see in the estimate of the constant C in the regression 



(4.6 ), we get a positive contribution to the risk premium. There will also be a term involving 
time to maturity, which will converge to zero in the long end of the futures curve. This part 
of the risk premium may contribute both positively or negatively. The CARMA factor is 
thus giving a positive risk premium for contracts, which start delivering reasonably soon. 
Since this factor is accounting for the short term variations, and in particular the spike risk 
of the spot, we may view this as a result of consumers and retailers hedging their price risk 
and, therefore, accepting to pay a premium for this. This conclusion is in line with the 
theoretical considerations of Benth, Cartea and Kiesel 0], who showed - using the certainty 
equivalence principle - that the presence of jumps in the spot price dynamics will lead to a 
positive risk premium in the short end of the futures curve. Bessembinder and Lemmon [8] 
explain the existence of a positive premium in the short end of the futures market by an 
equilibrium model, where the skewness in spot prices induced by spikes is a crucial driver. 



As we know from the third summand in Equation (4.17), the risk premium is dominated by 



a linear trend for most times to maturity except very short ones. In the latter case, the first 



summand of Equation (4.17) is dominating, leading to a small exponential decay when time 
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risk premium (base) risk premium (peak) 
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time to maturity time to maturity 



Figure 10: The estimated risk premium R pr (red) and the empirical risk premium R pr (black), 
for base load data (left) and peak load data (right). 



to maturity tends to 0. A plot of the empirical risk premium versus the theoretical one is 



given in Figure 10 



We see that for the base load contracts the positive risk premium in the short end of the curve 
is not that pronounced, however, it is detectable. The risk premium is negative for contracts 
starting to deliver in about two months. On the other hand, the peak load contracts have a 
clear positive risk premium, which changes to a negative one for contracts starting to deliver 
in about four months. This form of the risk premium is in line with the analysis of Geman 
and Vasicek [22] • Interesting here is the difference between base and peak load contracts. 
Base load futures have a longer delivery period than peak loads, since they are settled against 
more hours. This means that extreme prices are more smoothed out for base load contracts. 
The sensitivity towards spikes are even more pronounced in peak load contracts, since they 
concentrate their settlement for the hours, where typically the extreme spikes occur, and 
ignore to night hours where prices are usually lower and thus would smooth out the spikes. 
Hence, peak load contracts are much more spike sensitive than base load contracts, which we 
see reflected in the risk premium having a bigger and more visible positive part in the short 
end of the futures curve. The study of Longstaff and Wang (27] on the PJM (Pennsylvania, 
New Jersey and Maryland) market shows that the risk premium may vary over time, and 
indeed change sign. Their analysis is performed on hourly prices in the balancing market 
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as being the spot, and the day ahead hourly prices as being futures contracts. Hence, the 
analysis by Longstaff and Wang [27] is valid for the very short end of the futures curve. 

Our findings for the EEX are in line with the empirical studies in Benth et al. [I], which 
applies a two factor model to analyse the risk premium in the EEX market. Their model 
consists of two stationary processes, one for the short term variations, and another for 
stationary variations mean-reverting at a slower speed. Their studies confirm a change in 
sign of the risk premium as we observe for our model. Moreover, going back to Lucia and 
Schwartz |28J . they find a positive contribution to the risk premium from their short-term 
variation factor, when applying their analysis to NordPool data. This shows that also in 
this market there is a tendency towards hedging of spike risk in the short end of the futures 
curve. On the other hand, our results for the EEX market are at stake with the findings 
in Kolos and Ronn [26J. They perform an empirical study of many power markets, where 
they estimate market prices of risk for a two-factor Schwartz and Smith model. In the EEX 
market, they find that both the short and the long term factors contribute negatively]^] for 
the case of the EEX market. However, interestingly, the PJM market in the US, which is 
known to have huge price variations with many spikes observed, they find results similar 
to ours. We find our results natural in view of the spike risk fully accounted for in the 
short term factor, and the natural explanation of the hedging pressure from producers in the 
non-stationary factor. Our statistical analysis also strongly suggest non-Gaussian models for 
both factors, which is very different to the Gaussian specification of the dynamics in Kolos 
and Ronn [26]. 



6 Conclusion 

In this paper we suggest a two-factor arithmetic spot model to analyse power futures prices. 
After removal of seasonality, a non-stationary long term factor is modelled as a Levy process, 
while the short term variations in the spot price is assumed to follow a stationary stable 
5 In their paper, the signs are positive due to the choice of parametrization of the market price of risk 
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CARMA process. An empirical analysis of spot price data from the German power exchange 
EEX shows that a stable CARMA processes is able to capture the extreme behavior of 
electricity spot prices, as well as the more normal variations when the market is in a quite 
period. 

As in Lucia and Schwartz |28j we use a combination of a deterministic function and a non- 
stationary term to model the low frequency long term dynamics of the spot. Empirical data 
suggests that futures curves and spot prices are driven by a common stochastic trend, and 
it turns out that this is very well described by a normal inverse Gaussian Levy process. 
This leads to realistic predictions of the futures prices. Moreover, a CARMA(2,1) process is 
statistically the best model for the short term variations in the spot dynamics. 

We apply the Esscher transform to produce a parametric class of market prices of risk for the 
non-stationary term. The a-stable Levy process driving the CARMA-factor is transformed 
into a tempered stable process in the risk neutral setting. The spot price dynamics and the 
chosen class of risk neutral probabilities allow for analytic pricing of the futures. A crucial 
insight in the futures price dynamics is that the stationary CARMA effect from the spot 
price is vanishing for contracts far from delivery, where prices essentially behave like the 
non-stationary long-term factor. 

We propose a statistical method to calibrate the suggested spot and futures model to real 
data. The calibration is done using spot and futures data together, where we applied futures 
prices in the far end of the market to filter out the non-stationary factor in the spot. We 
choose a threshold for what is sufficiently "far out" on the futures curve by minimizing the 
error in matching the theoretical risk premium to the empirical. In this minimization over 
thresholds, we need to re-estimate the whole model until the minimum is attained. Since 
a-stable processes are not in L 2 , we introduce a robust L 1 -filter in order to recover the states 
of the CARMA process required for the estimation of the risk premium. 

Our model and calibration technique is used on spot and futures data collected at the EEX. 
Moreover, in order to gain full insight into the risk premium structure in this market, we 
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study both peak load and base load futures contracts with delivery over one month. The 
base load futures are settled against the hourly spot price over the whole delivery period, 
while the peak load contracts only deliver against the spot price in the peak hours from 8 
a.m. to 8 p.m. on working days. Our model and estimation technique seem to work well in 
both situations. 

We find that the base load futures contracts have a risk premium which is close to linearly 
decaying with time to delivery. The risk premium is essentially governed by the long term 
factor. There is evidence of a positive premium in the short end of the futures curve. For 
peak load contracts, which are much more sensitive to spikes, the positive premium in the 
short end is far most distinct, but also here the premium decays close to linearly in the long 
end of the market. These observations are in line other theoretical and empricial studies 
of risk premia in electricity markets, which argue that the risk premia in power markets 
are driven by hedging needs. Our findings also show that we should have an exponential 
dampening of the premium towards maturity, resulting from the CARMA factor of the spot. 
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